nanopyx.core.transform.sr_error_map

  1# adaptation of https://github.com/HenriquesLab/NanoJ-eSRRF/blob/master/src/nanoj/liveSRRF/ErrorMapLiveSRRF.java
  2
  3import numpy as np
  4from scipy.ndimage import gaussian_filter
  5from scipy.optimize import (
  6    brent,  # https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.brent.html
  7)
  8from scipy.stats import (
  9    linregress,  # https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.linregress.html
 10)
 11from skimage.transform import resize
 12
 13from ..analysis.pearson_correlation import pearson_correlation
 14
 15
 16class ErrorMap:
 17    def __init__(self):
 18        self._vRSE: float = 0
 19        self._vRSP: float = 0
 20        self._alpha: float = 0
 21        self._beta: float = 0
 22        self._sigma: float = 0
 23
 24        self.im_ref: np.ndarray = None
 25        self.im_sr: np.ndarray = None
 26        self.im_sr_intensity_scaled_blurred: np.ndarray = None
 27        self.imRSE: np.ndarray = None
 28
 29    def optimise(self, imRef: np.ndarray, imSR: np.ndarray, fixedSigma=0) -> None:
 30        self.im_ref = imRef
 31        self.im_sr = imSR
 32
 33        magnification = imSR.shape[0] / imRef.shape[0]
 34        assert magnification == imSR.shape[1] / imRef.shape[1]
 35
 36        if magnification > 1:
 37            imRef = resize(imRef, imSR.shape, order=3, preserve_range=True)
 38
 39        self.imRefMagnified = imRef
 40
 41        max_sigma_boundary = (
 42            4 / 2.35482
 43        ) * magnification  # this assumes Nyquist sampling in the ref image
 44
 45        sigma_linear = fixedSigma * magnification
 46        if fixedSigma == 0:
 47            sigma_linear = brent(
 48                sigma_function_to_optimize,
 49                args=(imRef, imSR),
 50                brack=(0, max_sigma_boundary),
 51                maxiter=1000,
 52            )
 53
 54        if abs(sigma_linear - max_sigma_boundary) < 0.0001:
 55            print("RSF constrained, as no good minimum found")
 56
 57        # GET ALPHA AND BETA
 58        alpha, beta = calculate_alpha_beta(sigma_linear, imRef, imSR)
 59        self._alpha = alpha
 60        self._beta = beta
 61        self._sigma = sigma_linear
 62        self.im_sr_intensity_scaled_blurred = gaussian_filter(
 63            imSR * self._alpha + self._beta, self._sigma
 64        )
 65        self.imRSE = np.abs(self.im_sr_intensity_scaled_blurred - imRef)
 66        self._vRSE = np.mean((self.im_sr_intensity_scaled_blurred - imRef) ** 2) ** 0.5
 67        self._vRSP = pearson_correlation(self.im_sr_intensity_scaled_blurred, imRef)
 68
 69    def getRSE(self) -> float:
 70        return self._vRSE
 71
 72    def getRSP(self) -> float:
 73        return self._vRSP
 74
 75    def get_sigma(self) -> float:
 76        return self._sigma
 77
 78
 79def calculate_alpha_beta(sigma: float, imRef: np.ndarray, imSR: np.ndarray) -> tuple:
 80    """Gaussian blurs imSR image and calculates linear regressino again imRef
 81
 82    Args:
 83        sigma (float): gaussian blur sigma
 84        imRef (np.ndarray): reference image (generally a difraction limited equivalent)
 85        imSR (np.ndarray): super-resolution image
 86
 87    Returns:
 88        tuple[float, float]: alpha and beta for linear regression
 89    """
 90    imSRBlurred = gaussian_filter(imSR, sigma)
 91    slope, intercept, r, p, se = linregress(imSRBlurred.ravel(), imRef.ravel())
 92    return slope, intercept
 93
 94
 95def sigma_function_to_optimize(
 96    sigma: float, imRef: np.ndarray, imSR: np.ndarray
 97) -> float:
 98    alpha, beta = calculate_alpha_beta(sigma, imRef, imSR)
 99    im_sr_intensity_scaled_blurred = gaussian_filter(imSR * alpha + beta, sigma)
100    rmse = np.mean((im_sr_intensity_scaled_blurred - imRef) ** 2) ** 0.5
101    return rmse